arXiv:1507.07373v2 [physics.comp-ph] 28 Jul 2015 


ATLAS: A Real-Space Finite-Difference Implementation of Orbital-Free Density 

Functional Theory 


Wenhui Xuecheng Shao®’', Chuanxun Su®, Yuanyuan Zhou®, Shoutao Zhang®, Quan Li®, Hui Wang®, Lijun 
Zhang’’, Maosheng Miao‘’, Yanchao Wang®’*, Yanming Ma®’* 

^State Key Laboratory of Superhard Materials, Jilin University, Changchun 130012, China 
^College of Materials Science and Engineering, Jilin University, Changchun 130012,China 
^Department of Chemistry and Biochemistry, California State University Northridge, 18111 Nordhojf St., Northridge, CA 91330, USA 


Abstract 

Orbital-free density functional theory (OF-DFT) is a promising method for large-scale quantum mechanics simula¬ 
tion as it provides a good balance of accuracy and computational cost. Its applicability to large-scale simulations 
has been aided by progress in constructing kinetic energy functionals and local pseudopotentials. However, the 
widespread adoption of OF-DFT requires further improvement in its efficiency and robustly implemented software. 
Here we develop a real-space finite-difference method for the numerical solution of OF-DFT in periodic systems. 
Instead of the traditional self-consistent method, a powerful scheme for energy minimization is introduced to solve 
the Euler-Lagrange equation. Our approach engages both the real-space finite-difference method and a direct energy- 
minimization scheme for the OF-DFT calculations. The method is coded into the ATLAS software package and 
benchmarked using periodic systems of solid Mg, Al, and AlsMg. The test results show that our implementation can 
achieve high accuracy, efficiency, and numerical stability for large-scale simulations. 

Keywords: Orbital-Free Density Functional Theory; Quantum mechanical; Real-Space Representations; local 
pseudopotentials 


1. Introduction 

Computational simulation is a powerful tool for 
predicting material properties and understanding the 
physics underlying experimental observations. |[l|] Reli¬ 
able simulation relies on advanced computational the¬ 
ories and methods, and in recent decades many effi¬ 
cient approaches with different levels of accuracy have 
emerged to receive remarkable success; e.g., quantum- 
mechanical [QSS] and empirical potential methods. Jsl 

i 

Quantum mechanical approaches based on the Kohn- 
Sham (KS) density functional theory (DFT) 10, [3] allow 
accurate descriptions of materials’ properties, but are 
computationally demanding. They require evaluation 
of the kinetic energy term related to the computation of 
single-electron wave-functions. The calculation of elec¬ 
tron density needs to consider 3Ae degrees of freedom. 
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and thus the computational cost scales as a cubic rela¬ 
tion of 0{NI), 10] where Ne is the total number of elec¬ 
trons of the system. The large computational cost limits 
KS-DFT to only small systems with unit cells of up to 
only hundreds or thousands of atoms, ||3] which makes 
its traditional implementation unsuitable for studying 
complex systems (e.g., surfaces, interfaces, nanomate¬ 
rials, and biomaterials) or macro-scale features (e.g., 
grain boundaries, dislocations, and cracks) that require 
large-scale simulations using cells of tens of thousands 
or millions of atoms. JH [3] 

Over the past two decades, many linear scaling tech¬ 
niques have been developed in an effort to reduce the cu¬ 
bic scaling of traditional KS-DFT. iSi However, th^ 
depend on both the “nearsightedness” principled |9|] 
and the concept of “locality”, m and therefore scale 
linearly only for systems containing a large number of 
atoms. There is an unavoidable crossover between cu¬ 
bic and linear scaling. inii Moreover, linear scal¬ 
ing requires a band-gap structure or localized electronic 
structure, and appears not to function for metals. mi 

Alternative approaches related to parameter fitting 
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have therefore been designed using empirical inter¬ 
atomic potentials. fin These simulations require much 
less computational cost and are computationally ca¬ 
pable of dealing with macro-scale problems, but they 
suffer notable shortcomings in terms of accuracy and 
transferability, ill More severely, these simulations 
completely neglect the properties of electrons, which 
are fundamentally important to various aspects of chem¬ 
istry and physics. There is therefore an urgent need for a 
reliable quantum-mechanics-based method able to per¬ 
form large-scale simulations. 

Orbital-free (OF) DFT,!!! [3 [3 Q [3 is poten¬ 
tially an efficient theory for large-scale quantum me¬ 
chanical simulations. The total energy within the OF- 
DFT scheme is expressed as an explicit functional of 
electron density in the Hohenberg-Kohn theorem, 13] 
and there is no need to deal with wave functions. Here 
the electron density, a simple function with three de¬ 
grees of freedom, can uniquely determine the ground- 
state properties of a many-electron system. As such, the 
computational cost of OF-DFT scales quasilinearly with 
the number of atoms in the system, providing substan¬ 
tial advantages in numerical simplicity and efficiency 
for large-scale simulations. IlllflzflMfldl fl^ 

The drawbacks of OF-DFT include two challenges to 
its realistic treatment of the kinetic energy density func¬ 
tional (KEDF) and ion-electron interactions. Ilisl] First, 
the kinetic energy term is a sole functional of the elec¬ 
tron density function; its construction determines the 
accuracy of OF-DFT, and significant mogress has been 
made in the last two decades (Refs, [[[li 


19, 


21, 22, 23, 24, 25, 2^ 13]]) with the proposal of various 
encouraging KEDFs. These functionals have been suc¬ 
cessfully tested in many systems, including nearly free- 
electron-like systems and semi-conductors. 113113] 
Second, the lack of particular orbitals leads the ion- 
electron interaction to be described only by the lo¬ 
cal pseudopotentials (LPPs), and previous studies have 
sought to construct various LPPs. Empirical |[l3 [idl 
[13, Ell 112] and “bulk derived” LPPs lIssI] have been de¬ 
veloped and successfully applied to various metals and 
semiconductors, m Esh but a notable challenge is the 
construction of a good LPP from an existing non-LPP 
without appeal to any bulk- or aggregate-system KS 
calculations. fill] We recently developed an optimized 
effective potential (OEPP) scheme to construct full first- 
principles LPPs from existing non- LPPsi3i to enhance 
the transferability of the pseudopotential. Our OEPP 
worked well for a large number of elements, and the 
transferability of the LPP was found to be an intrinsic 
property of elements. 

The ground-state energy Emin and electron density p 


in OF-DFtOEI can be obtained by minimizing the 
total energy E\p\ of the system with respect to the trial 
electron density p. The following minimization equa¬ 
tion is non-linear and multidimensional: BSl] 


Emin = min < E[p] 


-p( I p(r)dr- 
Jq 


W,);p>0 , (1) 


where p is a Lagrange multiplier used to enforce the 
constraint that the total number of electrons Ne is con¬ 
served, and Q is the whole space for the simulation. 

Currently, there are two main procedures for solving 
Eq. O. Ilisl] The first is to seek a direct solution of the 
equation by minimizing the total energy with respect to 
the electron density p.JSfl [STI [35 


38h The well-known 

Ei 


way, 


39 


and 


PROFESS code was built this 
has been successfully used to investigate many large- 
scale problems. ||4T1 l42l l43n The other procedure is to 
transform Eq. o into an ordinary KS-like equation 
that can be solved self-consistently by any KS com- 
puter pro gram. iH liM 1^ However, recent research fl3 
47l ksh has shown that the iterative self-consistent pro¬ 
cedure for OF-DFT does not work properly for large 
systems. Moreover, the non-convergence problem is 
not solved, and the underlying reason for this remains 
unclear. lO 

In this work, a real-space finite-difference method for 
solving the OF-DFT Euler-Lagrange equation (Eq. O) 
for periodic systems is developed by direct minimiza¬ 
tion. As shown previously, a real-space finite-difference 
method provides three obvious advantages: EH 

E2? ] (i) the method is independent of any basis, simpli¬ 
fying its implementation: D30 (ii) a real-space method is 
advantageous for large-scale parallel calculations due to 
its avoidance of the fast Fourier transform (FFT) method 
for the reciprocal-space approach and the serious draw¬ 
back in the need for “all-to-all communication” ||3] dur¬ 
ing parallel calculation; and (iii) there is no barrier to 
switching between a periodic and a non-periodic sys¬ 
tem for a real-space approach. Particularly, we find that 
the finite-difference method is computationally more ef¬ 
ficient in dealing with the Laplace, gradient, and diver¬ 
gence operators than the FFT-based method. 

Our method is coded into Ab initio orbiTaL-free den¬ 
sity functionAl theory Software (ATLAS) and is bench- 
marked in periodic systems of Mg, Al, and AlsMg. Our 
current implementation of OF-DFT is shown to be nu¬ 
merically accurate, stable, and efficient. 

The remainder of this paper is organized as follows. 
Section II gives the theory. The OF-DFT differential 
equation is presented for illustration, followed by de¬ 
tailed real-space representations of the finite-difference 
method and the direct energy minimization method to 
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obtain the ground-state electron density. Section III re¬ 
ports testing results on Mg, Al, and AlsMg crystals to 
demonstrate the computational accuracy, efficiency, and 
stability of the procedure. Finally, conclusions are pre¬ 
sented in Section IV. 


2. Theory and Background 


2.7. OF-DFT Theory 

In the OF-DFT scheme, fl^fl^ the ground state to¬ 
tal energy of an A-electron system in a local external 
potential Vextir) is a functional of electron density p(r): 


/ 


p{r)dr = Ne. 


( 2 ) 


The total energy functional E[p\ (atomic units, a.u., 
are used throughout the paper) can be written as 
follows: iilQ 


E^^\p\ =T\p\ + Eh\p\ + Exc\p\ + 


f 


Vext(r)p(r)dr 




(3) 


The third non-local term has the form 

E„i[p] = J'p‘‘(r)k[p(r),p(r),r,r']f/(r')drdr'. (7) 

This non-local term Tnilp] is Taylor expanded to achieve 
quasilinear scaling with system size via FFT. ll^l^ 
Our approach decomposes the KEDFs into the ATy^ 
term and a remaining term Tglp] (which, when A = 1, is 
the Pauli term lfl^ l: 


T[p]=ATyw{p] + Te{pl ( 8 ) 


The vW kinetic potential term can be evaluated as fol¬ 
lows: 


STywlp] 1 


(9) 


The finite-difference method is adopted here to obtain 
the vW kinetic potential term instead of FFT, as imple¬ 
mented in PROFESS. !^ 

Instead of minimizing directly over the electron 
density, we rewrite the total energy as the functional of 
0 = ^Jp. Taking the constraint of Eq. Q into account 
by Lagrange’s multiplier method, we define 


where Eulp] is the Hartree electron-electron repul¬ 
sion energy, Exc[p\ the exchange-correlation energy, 
and Vextir) the external potential representing the ion- 
electron interaction as given by local pseudopotentials 
in the OF-DFT scheme. Ei-i is the interaction en¬ 
ergy between ions, which is dealt with in our imple¬ 
mentation using Ewald summation. 10 The ki¬ 

netic energy of the non-interacting electrons {T[p\) is 
an explicit functional of electron density. These KEDFs 
can be roughly categorized into two general types: lo¬ 
cal/semilocal and nonlocal. The former naturally scale 
linearly with system size, while the later scales quadrat- 
ically owing to its double integral. Our program im¬ 
plements both local/semilocal KEDFs 101^ an d the 
nonlocal Wang-Govind-Carter (WGC) KEDfS. The 
WGC KEDF can be written as follows: 


nWGC 


[p] = Ttf{p\ + Tywlp] + TnAp], 


(4) 


where Ttf{p\ is the Thomas-Fermi term assuming the 
limit of a uniform electron gas. It takes the form 

Ttf{p\ = Ctf f p(r)dr, (5) 

Jq 

where Ctf = Tvwlp] is the von Weizsacker 

(vW) term designed for a single-orbital system: 

Fwlp] = £ FT) V^j FT)dr. (6) 


L{<p] = E' 


OFv 


X 


p \ I (p {r)dr - N 


then the gradient of L with respect to 0 is 
dL[0] SE^^[p] dp 


2p(f> 


where 


d0 Sp 6(p 

= 2{H<p-p<p}, 


H,p = --VFr) + Veff(r)m- 


( 10 ) 


This equality is derived from Eq. da and 

^ ^ _ dTelp] dEulp] SExclp] 

^effVn - —7 - + 7 - + 7 — 

dp dp dp 

= Voir) + VH{r) + yxc(r) + V.„(r). 

The variational principle requires dL/d0 
leads to the Euler-Lagrange equation: 


( 11 ) 

( 12 ) 

+ V,o„(r) 

(13) 
0, which 


H(p = pcp. (14) 

This is a Schrodinger-like equation, EiiliiSQ 
but much simplified as there is only one “orbital” with 
the minimal eigenvalue. This equation can then be re¬ 
solved by minimizing the OF-DFT total energy with re¬ 
spect to ^|p. 


ilESfis 
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2.2. Real-Space Representations 


Real-space calculations are performed on grids, in 
which the values of the electron density distribution and 
effective potential are given on discrete Cartesian grid 
points. The real-space finite-difference expansion trans¬ 
forms the kinetic energy operator into a spare matrix, 
which has nonzero elements only in the vicinity of the 
leading diagonal. Isil [12] The general form of the Lapla- 
cian with a Cartesian grid can be expressed as follows: 


N 

V'^(p(xi,yj,Zk) = ^ Cn^{xi + nhx,yj,Zk) 

n=-N 

N 

^ Cn(p(xi,yj-\-nhy,Zk) 

n=-N 

N 

^ Cn(/>(xi,yj,Zk + nhz), (15) 

n=-N 


where N is the order of the finite difference expansion; 
hx, hy, and are the grid spacings in the v, y, and z di¬ 
rections, respectively; and the coefficients are avail¬ 
able in Ref. lUll[fill]. However, the Cartesian grid is in¬ 
compatible with the periodicity of a non-orthorhombic 
unit cell. A new high-order finite-difference method for 
a non-orthorhombic grid has been proposed and suc¬ 
cessfully applied to periodic systems. 11621] We adopt it 
here. The general form of the Laplacian operator for a 
non-orthorhombic grid is as follows: 





(16) 


We represent the Laplacian by a combination of deriva¬ 
tives along six nearest-neighboring v, directions: three 
original at (i= 1,2,3) directions and three additional 
nearest-neighboring directions, where at are the lattice 
vectors in real space. For the coefficient, refer to Ref. 

Note that the H matrix (Eq. (O) is a spare matrix 
whose nonzero elements are confined within a diago¬ 
nal band, and the extent of the nonzero elements in off- 
diagonal positions depends on the order of the finite dif¬ 
ference expansion. The Hartree potential is determined 
by solving the Poisson equation: 

V"VH(r) = -47r|p(r)-po(r)], (17) 


where po(r) is the average electron density of the sys¬ 
tem. For infinite periodic systems, we encountered the 
divergent problems on the ion-electron, ion-ion, and 
electron-electron interaction energies arising from the 
long-range Coulomb interaction -Z/r. Fortunately, the 



r (a.u.) 

(a) 



r (a.u.) 
(b) 


Figure 1: Local optimized effective pseudopotential of Mg (a) and A1 
(b) in real space. 
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divergent problem converts into the singularity problem 
at g = 0 in reciprocal space. For a charge-neutral sys¬ 
tem, the singularity for Hartree electron-electron poten¬ 
tial can be exactly canceled by adding up singularities 
encountered in the electron-ion and ion-ion potentials, 
and can therefore be neglected. |[^[35|] The Hartree po¬ 
tential Vuir) can thus be obtained as follows: 


VH{r) = FFT\^p(G)) 

|CrF 


(p(G = 0) = 0), (18) 


where p(G) is the electron density in reciprocal space, 
and FFT is a reverse FFT transform. The ionic term 
Vioni't') in Eq. (O can be constructed from Viodr) (i.e., 
LPPs). We use our developed OEPP for LPPs. OEPPs 
for both Mg and A1 are shown in Eig. |l(a)| and |l(b)[ 
respectively. The theory of constructing OEPP is pre¬ 
sented elsewhere. S Eor a periodic system, the ionic 
potential receives contributions from an infinite number 
of atoms, leading to a divergent summation of the long- 
range Coulomb term. To seek a solution, as mentioned 
above, the pseudopotential Vionjr) is then expressed in 
reciprocal space as follows: ( tI 

VioniG) ^ ^ r 


yioniT)exp{iG • r)dr 


1 




'cell 


•we 

Y,s\g)vi^{G), 

K=\ 


(19) 


where Utype is the number of atomic species, and for 
each atomic species k there are identical atoms at 
positions = 1,^^, and Tlceii is the unit cell vol¬ 

ume. The structure factor S (G) for each atomic species 
Khl 


S\G) = V expiiG • 


r- 


( 20 ) 


and the form factor ViociG) is 

KciG) = f yi,{r)exp{iG • r)dr. (21) 

J allspace 

The spherical symmetry of LPP allows the 3D Eourier 
transform to b^represented by a ID radial Eourier 
transform:! 


KM) 


0 


N 2^in{gr) j 

KM)^ - dr 

LUC 


= C(^)-— (^^0), (22) 

where the non-Coulomb potential V^cig) in recipro¬ 
cal space can be written as 

Vi(*) = 4. (23) 

Jo ^ 8^ 


At g = 0, the Coulomb interaction is canceled as de¬ 
scribed above, leading to 


rKut 

VL(g = 0) = 47r Mv^Jr) + Zr)dr, 
Jo 


(24) 


where rcut is the cutoff of core radii. In our implementa¬ 
tion, Vioc(^) is equal to -Z/r when r > Vcut- 

Eor a given grid spacing h, the size of the grid points 
can be determined as where 


N, = 


(25) 


Eor a given structure, the wave vector G is determined 


by 


G{ni,n2, ^ 3 ) = nibi ^ 2^2 + ^ 3 ^ 3 , 


(26) 


where bt (i = 1,2, 3) are the primitive vectors in recip¬ 
rocal space, and = (0,1,2 • • -Nt) are integers. 

Einally, the real-space local ion-electron pseudopo¬ 
tential Vioni^) can be calculated by an EET: 


V.-o„(r) = FFT(Vton(G)). 


(27) 


As mentioned above, all the physical quantities in the 
real-space finite-difference formalism can be directly 
represented on discretized grid points with a uniform 
interval. 11621] Einally, Eq. STSi can be expressed as the 
following discretized expression: 

H(p(ui,Vj,Wk) = 

A ^ 

--[ ^ Cn,(l){uiM^nniKvjM^i2niKwk + ^unih) 
ni=-N 
N 

+ ^ Cn2(p(Ui-\-^2in2KVj-\-^22n2KWk-\-^23n2h) 
n2=-N 
N 

+ Z Cn,<t>i^i + ^3in3h,Vj+^32n3h,Wk+^33n3h) 
n,=-N 
N 

+ E Cn^ct){Ui + ^4in4h,Vj + ^42n4h, Wk + ^43«4^) 

n/^=-N 

N 

+ ^ Cn^(p{UiM ^sinsKVj M 

n5=-N 

N 

ne=-N 

HVion(Ui,Vj, Wk) + VH(Ui,Vj, Wk) + Vxc(l^hVj, m) 

MVeiui^Vj, Wk)] X (f)(ui,Vj, Wk), (28) 

where C„. = fiCm, N is the order of the finite-difference 
expansion. The choice of ^ as -1, 0, or 1 depends on the 
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lattice vectors of the crystal. Given that the Laplacian 
operator extends to only a few neighbors around each 
grid point, Eq. (1^ is a sparse matrix. We solve it here 
to obtain the minimum energy [p] by implementing 
an energy-minimization scheme. Previous works have 
shown that the Truncated Newton (TN) method ic is 
one of the most efficient; [ssl therefore, it is 
employed in our approach. 

2.3. Algorithm for Energy Minimization 

Our scheme selects ^fp as the variable to minimize 
the total energy. The flow chart of ATLAS code is 
shown in Fig. O The procedure followed comprises 
three major steps. First, an initial guess of the electron 
density po is required for a trial solution of (p^. Note that 
the initially guessed electron density is derived from 
the model of a homogeneous electron gas (the average 
charge density of the system). 

Second, the ground-state electron density is obtained 
by minimizing the total energy using the TN method 
based on the initial guess. The TN algorithm for energy 
minimization consists of two iterations: an outer iter¬ 
ation that approximates the descent direction p as the 
direction for minimizing energy, and an inner iteration 



that determines the s^ size ^ by a line search to ensure 
an energy decrease. 0, [ssl [s^] The details of this 
step include three procedures. 

(i) According to the TN scheme, the search direction 
\pk) at iteration k is simply determined by the quantities 
of the current iteration. \pk) can be written as follows: 


\Pk) = -AZ\gk), 


(29) 


where \gk) is the gradient of L^(0) according to Eq. ([TT]) . 
which can be written as 


with 


\gk) = 2 {Hk\(pk} - Eklfk )}, 

{(pk\Hk\(pk) 


Ek = 




Ak is the approximate Hessian matrix of Lk'. 
. 8^Lt 


d0(r)d0(r') 

As in previous works, [3^ 35] we rewrite Eq. 
linear equation to determine \pky. 

M\Pk) = 


(30) 


(31) 


(32) 


as a 


(33) 


This equation can be solved using the linear conju¬ 
gate gradient method. We compute Ak\p) using the 
first-order finite-difference approximation rather than 
attempt the explicit evaluation of A^: 


M\p) 


\g{f + ep)) - \g{f)) 


(34) 


This ensures that the computational cost of our approach 
are linear scaling. 

(ii) The step size 6k is determined by line search with 
the normalization constraint of \(pk+i}- IPk) is further or- 
thogonalized to \(pk} and normalized to Ng. 


\(Pk > = \Pk} -— IPk) 


Ne 




Ky 


(35) 


(36) 


\(l>k+i) then is updated by 

l<^i+i> = \<Pk) cos(0k) + \<pf) sin(6»i), (37) 

where the value of 0ic is determined by line search!^ 


67L l68[] with the Wolfe conditions to ensure it lies toward 


lower energy. in 
Ok 


minE [fk+iix, 0)] 

6 


(38) 
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(iii) When the step size 6k is determined, the new 
electron density can be derived by pk+i = (pl+y The 
process is repeated until both the gradient of Lagrange 
\g) and the variation of total energy are smaller than the 
given tolerances. 

Finally, the third step involves calculating the total 
energy or other related physical quantities of a given 
structure from the ground-state electron density. 


3. numerical results 

We consider here three bulk systems of Mg, Al, and 
AlsMg to benchmark the above formalism (as imple¬ 
mented in our ATLAS code) for accuracy and compu¬ 
tational efficiency. The calculation employs the local 
density approximation (LDA) for electron exchange and 
correlation as parametrized by Perdew and Zunger.ll^ 
S The local pseudopotentials of Mg and Al are con¬ 
structed by our OEPP scheme for their respective elec¬ 
tronic configurations of 3s^3p^ and 3s^3p^. The core 
cutoff radii are 2.6 a.u. for Mg and 2.2 a.u. for Al. 

3.1. Tests of Real-Space OF-DFT Convergence 

Our real-space finite-difference implementation of 
OF-DFT has two controllable parameters that critically 
influence the accuracy of the calculations: the order of 
finite difference expansion and the grid spacing h. These 
parameters are chosen depending on the convergence 
test of the total energies of the systems. The grid spac¬ 
ing in real space is related to the plane-wave cutoff en¬ 
ergy {Ecut = 71^12hf) in reciprocal space. Here we show 
in a real application how to choose the values of these 
parameters. We run ATLAS code on calculations of to¬ 
tal energy for bulk Mg with a body-centered cubic (bcc) 
lattice. Fig. [3] shows that a fourth-order finite-difference 
expansion and a grid spacing of 0.18 A are sufficient for 
a well-converged total energy (0.1 meV/atom). Similar 
results are also found for bulk Al. Therefore, these two 
values are adopted for all the following calculations on 
systems of bulk Mg, Al, and AlsMg. 

3.2. Computational Accuracy 

For benchmarking our ATLAS program, bulk prop¬ 
erties of hexagonal close-packed (hep) Mg and face- 
centered cubic (fee) Al and AlsMg are calculated 
and compared with those calculated by the CASTEP 
codej^ within KS-DFT. Our ATLAS calculations em¬ 
ploy the WGC formula for KEDF ll^ (with the fol¬ 
lowing parameters: y = 2.7, a = (5 -r V5)/6, and 
J3 = (5 - V5)/6) and the OEPP local pseudopoten¬ 
tial. Note that the WGC form of KEDF is known to 




(b) 


Figure 3: Effect of (a) grid spacing and (b) order of finite-difference 
approximation on the total energy of bulk Mg with a bcc lattice 


describe accurately various bulk properties of Mg and 
Al. For meaningful comparison, the CASTEP calcu¬ 
lations adopted the same OEPP local pseudopotentials 
as used in our method. The calculated equilibrium vol¬ 
umes, total energies, and bulk moduli are listed in Table 
[U Our ATLAS results are in excellent agreement with 
those obtained by CASTEP; the noticeable small dif¬ 
ferences stem from the difference between the kinetic 
energy terms used in the two codes. 

We now focus on the fundamental quantity of elec¬ 
tron density as calculated by the ATLAS and CASTEP 
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(C) (d) 

Figure 4: Contour plots of electron density calculated by ATLAS (solid line) and CASTEP (dotted line), (a) The (001) and (b) (Oil) planes of Al; 
(c) the (OOOl)and (d) (0110) planes of Mg. Electron density and lattice vectors are given in A and (a.u./lOO), respectively 


Table 1: Bulk properties obtained by OE-DFT and KSDET methods: 
equilibrium volume (Vq per atom in A^), total energy (Eq in eV per 
atom), and bulk moduli (Bq in GPa) 


systems 

methods 

Vo 

Eo 

Bo 

Mg 

KS 

22.023 

-24.588 

36.5 


OF 

22.225 

-24.577 

35.0 

Al 

KS 

18.029 

-56.799 

69.4 


OF 

18.435 

-56.801 

67.4 

AlsMg 

KS 

19.031 

-48.767 

55.2 


OF 

19.019 

-48.757 

57.5 


codes. Fig. |4]shows that the two calculations give essen¬ 
tially identical contour plots of electron density in the 
(001) and (011) planes for Al and the (0001) and (0110) 
planes for Mg, thus supporting the accurate implemen¬ 
tation of OF-DFT in ATLAS code. Further validation 
of our method is given by randomly generating ten dif¬ 


ferent structures of Mg using the CALYPSO software 
package ||^[^ as listed in Tabled and then calculat¬ 
ing their total energies using both codes with the same 
OEPP. The results (Fig. O show expected small energy 
differences between the two sets of data, but both cal¬ 
culations give essentially identical structure sequences 
in energy order, providing further confidence in the ro¬ 
bustness of our ATLAS code. To verify the correct¬ 
ness of our new implementation, we calculate the total 
energy of Mg, Al, AlsMg, and other more complex sys¬ 
tems (e.g., distorted structures of fee Mg with big cells 
containing atomic distortions following a frozen phonon 
at the smallest wave vectors and with the length of the 
longest cell vector varying from 14.2 to 142 A) using 
profess !^ with the same OEPP and KEDF (e.g.. 


TFTvW T = 1, 1/5, 1/9) as used in ATLAS. The results 
show that the energy difference obtained between AT¬ 
LAS and PROFESS is less than 0.1 meV/atom for all 
these systems. 
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Table 2: Details of ten random structures of Mg generated by CALYPSO 


Structure no. 

Space group (Number) 

-C- 

Lattice parameters (A) 

Wykoff position 

1 

P2i/c(14) 

a=2.9688 

4e -0.51477 0.74802 -0.24572 



b=4.7352 

4e -0.78771 0.52847 -0.60354 



c=lL5362 

112.3421 

2d 0.50000 -0.00000 -0.50000 

2 

P-6m2(187) 

a= b=5.6880 

lb 0.00000 0.00000 0.50000 



c=2.6768 

Id 0.33333 0.66667 0.50000 

3j 0.81547 0.18453 0.00000 

3 

Ima2 (46) 

a=4.2219 

4a 0.00000 0.00000 0.07444 



b=8.5195 

4b 0.25000 0.29658 0.01732 



c=8.3406 

4b 0.25000 -0.16647 0.26886 
4b 0.25000 0.55467 0.11780 

4b -0.25000 0.68084 0.26549 

4 

P3ml (156) 

a=b=5.6997 

la 0.00000 0.00000 0.23588 



c=5.3316 

la 0.00000 -0.00000 0.65630 
lb 0.33333 0.66667 0.52999 

Ic 0.66667 0.33333 0.65311 

3d 0.182460.817540.99311 

3d 0.50622 0.49377 0.25310 

5 

Pnn2 (34) 

a=4.6491 

4c 0.22362 0.16685 0.57664 



b= 11.0429 

4c 0.70660 0.63845 0.95432 



c=2.9217 

2a 0.00000 0.00000 0.08609 

6 

P4/mmm (123) 

a=5.9328 

Ic 0.50000 0.50000 -0.00000 



b=5.9328 

Id 0.50000 0.50000 0.50000 



c=4.2617 

41 0.32079 0.00000 -0.00000 
4m 0.29232 0.00000 0.50000 

7 

P4 (75) 

a=b=5.7234 

4d 0.70357 0.72197 0.85041 



c=4.5792 

la 0.00000 0.00000 0.33461 
la 0.00000 0.00000 0.85646 

2c 0.00000 0.50000 0.03534 

2c 0.00000 0.50000 0.49923 

8 

P222i (17) 

a=5.3893 

4e 0.69544 0.72370 1.34282 



b=4.7685 

2a 0.954510.00000 0.50000 



c=5.8367 

2a 0.41997 0.00000 0.50000 

2c 0.00000 0.67949 0.75000 

9 

P2 (3) 

a=6.2958 

2e -0.66003 -0.52140 1.47599 



b=7.6726 

2e -0.73000-0.98151 1.40517 



c=3.1136 

2e -0.15425 -0.13749 1.09614 



y8 = 94.1982 

2e -0.74070 -0.72248 0.96147 
lb -0.00000 -0.69511 0.50000 
Ic -0.50000 -0.34474 1.00000 

10 

Pmna (53) 

a=3.3393 

2b 0.00000 0.00000 0.50000 



b=7.5958 

4h 0.50000 0.30656 0.43343 



c=5.9138 

4h 0.50000 0.75004 1.12286 


3.3. Computational Efficiency 

Note that a prominent difference between our method 
and prior works is that the vW term is evaluated 
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Figure 5: Relative energy differences between the first structure and 
another nine structures of Mg generated by CALYPSO. The numbers 
along the horizontal axis correspond to the structures in Table|2] which 
lists their detailed structural information. 



0.0 1.0E7 2.0E7 3.0E7 4.0E7 5.0E7 6.0E7 7.0E7 

Number of grid points 


Figure 6: Timings (wall time) used to calculate vW kinetic potential 
in different numbers of grid points with both the finite-difference ex¬ 
pression and the FFT-based method. 


with a finite-difference expression instead of the FFT- 
based approach. The (wall) times for calculating vW 
kinetic potential for different sizes of grid points via 
fourth-order finite-difference expression and the FFT- 
based method using FFTwItI are shown in Fig. 1 
The finite-difference approach is clearly computation¬ 
ally more efficient than the FFT method, especially for 


Figure 7: Timings (wall time) using ATLAS to calculate the total en¬ 
ergy within the course of an electron-density optimization for systems 
of 10 to 10000 atoms in a simulated bcc Mg cell. The total time (blue 
line) is shown as the sum of the times for the ion-electron potential 
term (black line) and for all other potential terms and energy terms 
(green line). Also shown is the total time (wall time) cost for static 
energy calculation on systems of 2 to 240 atoms using CASTEP (red 
line) 


denser grid points. This is due to the different size de¬ 
pendence of the two methods. Assuming N is the num¬ 
ber of grid points, the computational cost of the finite- 
difference method is proportional to 0(N), whereas that 
of FFT is proportional to 0(NlogN). Note that our ap¬ 
proach shows a similar advantage in dealing with the 
generalized gradient approximation (GGA) kinetic po¬ 
tentials relating to the gradient and divergence opera¬ 
tors. 

A further test of the computational efficiency of the 
ATLAS software package is given in the analysis of bcc 
Mg. Fig. Ill shows the (wall) times for calculations of 
ion-electron potential terms, all other potential terms in 
Eq. (fTSl) . and the total energies within the course of an 
electron density optimization on systems containing 10 
to 10,000 atoms using a single processor. For compar¬ 
ison, single-processor calculations are also performed 
using the DFT code of CASTEP for systems containing 
up to 240 atoms. Both systems use the same exchange- 
correlation functional and OEPP. As expected, ATLAS 
shows a substantial advantage in computational effi¬ 
ciency over KS-DFT calculation. Note that the number 
of iterations to reach convergence (8-10) changes little 
with system size. Our method therefore shows strong 
potential applicability to large-scale simulation. 
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In fact, the computational cost of ATLAS shows 
quadratic scaling instead of linear, because the ion- 
electron potential term involving an explicit treatment 
of structure factors scales quadratically for a periodic 
system. 10 To avoid the quadratic scaling problem, a 
particle-mesh Ewald algorithm, which has linear scal¬ 
ing for the ion-electron term, will be implemented in 
ATLAS for periodic systems Note that for an 

isolated system there is no need to compute the structure 
factor. As a result, the ion-electron potential term will 
naturally show linear scaling. 

3.4. Numerical Stability 

The calculations for the Mg, Al, and AlsMg sys¬ 
tems show that ATLAS is numerically stable with 
TF/lvW and WGC KEDF. However, previous works ll^ 
in have indicated that procedures implementing most 
GGA KEDFs are numerically unstable. In the previous 
implementation, the numerical evaluation of gradient 
and divergence operators used the EFT. In our method, 
we adopted the finite-difference method to evaluate 
these operators for all the GGA KEDFs, [H] and 
tested our method using several systems (e.g., Al and 
Mg with fee, bcc, hep, and sc structures). The re¬ 
sults indicate that ATLAS with GGA KEDFs is also 
numerically unstable except for TFTvW, EOO, and P92, 
which is consistent with previous work. jUl There- 
fore, we believe that the numerical instabilities of most 
GGA KEDFs originate from unphysical electron den¬ 
sity produced by the singular and unphysical kinetic 
potential iQ [3 [13] during the process of optimizing 
the electron density. 

4. conclusion 

We developed an efficient ab initio method for the 
numerical solution of OF-DFT for large-scale simula¬ 
tions on periodic systems, and coded it into the AT¬ 
LAS software package. Our method employs the real- 
space finite-difference formulation and the scheme of 
energy minimization to yield both computational ac¬ 
curacy and efficiency for large-scale simulations. The 
performance of our method is well tested by designed 
static simulations for periodic systems of Mg, Al, and 
AlsMg, as well as comparison with data obtained by 
previous OF-DFT (PROFESS) and KS-DFT software 
packages (CASTEP). The results reveal that, except for 
the ion-electron term, the computational costs of the 
calculations of all other potential terms scale linearly 
with system size for periodic systems. Our future de¬ 
velopments of ATLAS code will focus on the imple¬ 


mentation of linear scaling particle mesh Ewald algo¬ 
rithms in an effort to achieve linear scaling on the ion- 
electron term.lir^. 13 0 more efficient algorithms for 
energy minimization, compatibility with non-periodic 
systems, parallel computing, and the evaluation of force 
and stress for ion and cell relaxations. We believe that 
ATLAS will become an alternative method for large- 
scale ab initio simulations. 
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